ON THE EXISTENCE OF THREE-DIMENSIONAL 
HYDROSTATIC AND MAGNETOSTATIC EQUILIBRIA 
OF SELF-GRAVITATING FLUID BODIES 



Daniele Galli 

INAF-Osservatorio Astrofisico di Arcetri 
Largo Enrico Fermi 5 
1-50125 Firenze, Italy 
galli@arcetri.astro.it 

ABSTRACT 

We develop an analytical spectral method to solve the equations of equilibrium for a 
self-gravitating, magnetized fluid body, under the only hypotheses that (a) the equation 
of state is isothermal, (6) the configuration is scale- free, and (c) the body is electrically 
neutral. All physical variables are represented as series of scalar and vector spherical 
harmonics of degree I and order m, and the equilibrium equations are reduced to a set of 
coupled quadratic algebraic equations for the expansion coefficients of the density and 
the magnetic vector potential. The method is general, and allows to recover previously 
known hydrostatic and magnetostatic solutions possessing axial symmetry. A linear per- 
turbation analysis of the equations in spectral form show that these basic axisymmetric 
states, considered as a continuos sequence with the relative amount of magnetic support 
as control parameter, have in general no neighboring nonaxisymmetric equilibria. This 
result lends credence to a conjecture originally made by H. Grad and extends early 
results obtained by E. Parker to the case of self-gravitating magnetized bodies. The 
only allowed bifurcations of this sequence of axisymmetric equilibria are represented 
by distortions with dipole-like angular dependence (I = 1) that can be continued into 
the nonlinear regime. These new configurations are either (i) azimuthally asymmetric 
(m = ±1) or (m) azimuthally symmetric but without reflection symmetry with respect 
to the equatorial plane (m = 0). It is likely that these configurations are not physically 
acceptable solutions for isolated systems, but represent instead the manifestation of a 
general gauge freedom of self-similar isothermal systems. To the extent that interstellar 
clouds can be represented as isolated magnetostatic equilibria, the results of this study 
suggest that the observed triaxial shapes of molecular cloud cores can be interpreted in 
terms of weakly damped Alfven oscillations about an equilibrium state. 
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1. Introduction 

Much of what we know of magnetostatic (MS) equilibria in plasma physics and astrophysics has 
come from studying systems with one ignorable coordinate: axisymmetric systems are one example. 
The symmetry associated to with one ignorable coordinate conveniently allows MS problems to be 
reduced to neat tractable forms. In the three-dimensional case, no one has ever been able to 
find nonaxisymmetric solutions of the MS equations for a self-gravitating fluid, and their actual 
existence has remained dubious (see Sect. 1.2). The goal of this paper is to explore the the existence 
of configurations of equilibrium of self-gravitating fluid bodies in the presence of a large-scale 
magnetic field without symmetry restrictions, with an application in mind to the study of the 
densest parts of molecular clouds (the so-called cloud cores), the sites of star formation. 

The motivation for this work is the following. The dense cores of interstellar molecular clouds 
are observed to be close to a state of virial equilibrium (balance of gravitational, thermal and 
magnetic energy (see e.g. Myers 1999), yet their shapes are distinctly non-spherical: the best fit 
to the projected axial ratio distribution is obtained with triaxial ellipsoids (Jones & Basu 2002; 
Goodwin, Ward-Thompson, &; Whitworth 2003). The question is then to consider the possible 
role of the anisotropic forces associated to the magnetic fields present in these cores (see e.g. 
Crutcher 2001) to allow the existence of shapes of non-trivial topology. 

The paper is organized as follows: in Sect. 2 we summarize the available results on the equi- 
librium of ideal (i.e. non viscous and non resistive) plasmas, discussing in particular the case of 
non-selfgravitating configurations; in Sect. 3 we formulate mathematically the problem and intro- 
duce non-dimensional self-similar variables; in Sect. 4 we describe the method of solution, based 
on the expansion of the angular part of the variables in scalar and vector spherical harmonics; in 
Sect. 5 and 6 we show that our method of solution recovers previously known axially symmetric 
hydro- and magnetostatic equilibria, respectively; in Sect. 7 we linearise the equations around these 
basic equilibria, and look for neutral modes (neighbouring equilibria) possessing a different symme- 
try than the basic states; finally, in Sect. 8 we summarize our results and discuss their implications 
for the observed shapes of molecular cloud cores. 



2. Ideal magnetostatic equilibria 

Consider a magnetized, isothermal, self-gravitating fluid body satisfying the ideal MS equations 

c s 2 Vp + P VV = ^(VxB) xB, (1) 

V • B = 0, (2) 

V 2 V = 4vrGp, (3) 

where B is the magnetic field, p is the density, V is the gravitational potential, G is the constant of 
gravitation and c s is the sound speed. With applications to interstellar molecular clouds in mind, 
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we have assumed here an isothermal equation of state. All known solutions of the set of equations 
(l)-(3) are characterized by a coordinate symmetry that reduces the number of variables from 
three to two or one. The symmetry associated with an ignorable coordinate allows MS problems 
to be reduced to a second-order elliptic partial differential equation (Dungey 1953), conventionally 
called the Grad-Shafranov equilibrium equation. These symmetric equilibria may be grouped into 
(i) axisymmetric (d/d<p = 0), (ii) cylindrically symmetric (d/dz = 0), and (Hi) helically symmetric 
systems (d/dip = kd/dz). As shown by Edenstrasser (1980), helical symmetry represents the 
most general admissible invariance property of the MS equations, with rotational and translational 
invariance as limiting cases. 

Applications of eq. (l)-(3) to the study of interstellar molecular clouds have largeley focused on 
axisymmetric magnetic configurations (e.g. Mouschovias 1976; Nakano 1979; Mestel & Ray 1985; 
Tomisaka, Ikeuchi,& Nakamura 1988; Barker & Mestel 1990). Cylindrically symmetric equilibria 
have also been extensively studied, originally in connection with the stability of galactic spiral arms 
(Chandrasekhar & Fermi 1953; Stodolkiewicz 1963), and later as models of filamentary clouds (Na- 
gasawa 1987; Nakamura, Hanawa, & Nakano 1995). In the latter context, magnetic configurations 
possessing helical symmetry have also been explored (Nakamura, Hanawa, & Nakano 1993, Fiege 
& Pudritz 2000a,b). An excellent introduction to the subject can be found in the monography by 
Mestel (1999). 



2.1. Parker's theorem 

For non-selfgravitating plasmas, the fundamental question as to the existence of MS equilibria 
of a more general symmetry, or with no symmetry at all (3-D MS equilibria), has been formulated 
several times (e.g. Low 1980, Degtyarev et al. 1985) but never properly answered. Grad (1967) 
conjectured that, with V = 0, only "highly symmetric" solutions of the system (l)-(2) should be 
expected, in order to balance the highly anisotropic Lorentz force with with pressure gradients and 
gravity, which are forces involving scalar potentials. In a fundamental paper, Parker (1972) proved 
rigorously the non existence of 3-D MS equilibria that are small perturbations of 2-D equilibria 
having translational symmetry. This result is conventionally referred to as Parker's theorem. 

On the basis of this result, Parker (1979) argued that realistic magnetic fields with no well 
defined symmetries must evolve in a genuinely time-dependent way, until all non-symmetric com- 
ponents of the field are destroyed by dissipation and reconnection and the topology becomes sym- 
metric, a process known as topological non- equilibrium of configurations lacking high degrees of 
symmetry (see also Tsinganos, Distler, &; Rosner 1984). The numerical simulations of Vainshtein 
et al. (2000) provide a striking illustration of this process. The problem with interstellar clouds is 
that ohmic dissipation times are of the order of ~ 10 15 yr, and therefore the kind of monotonic re- 
laxation to equilibrium envisaged by Parker (1979) can be ruled out. For a non-dissipative plasma, 
Moffatt (1985, 1986) has shown that stable non-axisymmetric magnetostatic equilibria of non-trivial 
topology do exist, but may require the presence of tangential field discontinuities (current sheets). 
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The mathematical method adopted in this paper, based on analytical functions and regular per- 
turbation expansions, cannot address the question of the existence of Moffatt-type equilibria in 
self-gravitating ideal plasmas. However, given the general character of Moffatt's conclusions, the 
question is considered again at the light of our results in Sect. 8. 

Does the presence of fluid motions modify this picture? Tsinganos (1982) found that Parker's 
theorem remains valid for steady dynamical (v 7^ 0) configurations possessing translational invari- 
ance, and stressed the analogy between this result and the familiar Taylor-Proudman theorem of 
hydrodynamics. However, Galli et al. (2001) found that a class of two-dimensional axisymmetric 
MHD equilibria (rotating, magnetized, self-gravitating singular isothermal disks) does have neigh- 
boring non-axisymmetric equilibrium states, provided the rotation speed becomes sufficiently high 
(supermagnetosonic). Thus, not only the presence of fluid motions, but also the geometry of the 
system, seems to play a crucial role in breaking the symmetry of the equilibrium state. 

As for the effect of gravity, Field (1982, quoted by Tsinganos et al. 1984) objected that the 
neglect of the constraining effect exerted by the plasma's self-gravity may severely limit the domain 
of existence of MS equilibria. In a series of papers (Low 1985; Bogdan & Low 1986; Low 1991), Low 
and collaborators have elaborated a general method to solve the equations of MS in the presence 
of an external gravitational field, but the problem presents considerable mathematical difficulties. 
In agreement with Grad's conjecture, Low concluded that applying some form of symmetry to 
the magnetic field is probably essential for the existence of equilibrium, in order to balance the 
highly anisotropic Lorentz force with pressure gradients and gravity, which are forces involving 
scalar potentials. In the special case of an imposed gravitational field, either uniform or due to a 
point mass, Low was able to find families of three-dimensional MS solutions. For a self-gravitating 
gas, Low (1991) showed that the problem can be reduced to the solution of two coupled partial 
differential equations for two unknown functions, but did not proceed further. 



3. The equations of the problem 

In the following we specialize to scale- free equilibria. Working in spherical coordinates (r, 6, ip), 
we assume that every physical quantity of the problem can be factorized in one function (power- 
law or logarithmic) of r times a function of 9 and <p only. This assumption allows a considerable 
simplification of the equations of the problem. We then show that the magnetic field must be 
poloidal in the Stratton-Chandrasekhar classification (Sect. 3.1), and we derive the governing set 
of non-dimensional equations (Sect. 3.2). Finally, we obtain the equations of the problem in the 
special case of axial symmetry (Sect. 3.3). 

The assumption of a power-law (or logarithmic) dependence on radius for the variables of the 
problem can be justified on observational and theoretical grounds. The observed density profiles 
of cloud cores approach in general a r~ 2 power-law behaviour in the outer parts, but are flatter 
near the cloud's centre (see e.g. Ward-Thompson, Motte, & Andre 1999). Idealised models that 
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assume a power- law behaviour of the density (and the intensity of the magnetic field, etc.) over 
all radii are justified by our expectation that the evolution of the observed cloud cores tends to 
a singular configuration. The driving process responsible for this evolution has been identified in 
the ambipolar diffusion of a weakly ionised gas, as originally proposed by Mestel & Spitzer (1956). 
Several numerical calculations of cloud evolution driven by ambipolar diffusion show that the density 
profile steepens to r~ 2 to arbitrary number of decades in radius as the singular state is approached 
(Fiedler & Mouschovias 1993; Ciolek & Mouschovias 1993, 1994; Basu & Mouschovias 1994, 1995). 
Thus, scale-free (singular) configurations provide realistic models for molecular cloud cores in the 
so-called "pivotal" state, i.e. on the verge of gravitational collapse. The price to pay for this 
simplifying assumption, is, of course, the introduction of an artificial singularity at the origin, 
where the density and the magnetic field diverge. 



3.1. Toroidal and poloidal fields 



Any solenoidal vector field B can be expressed as a linear combination of a certain basic 
toroidal (B t ) and poloidal (B p ) fields (Stratton 1941, Chandrasekhar 1961), given by 



r 



and 



V x 



V 



x r, 



x r 



(4) 



(5) 



where and ^ are arbitrary scalar functions of position. In addition to have zero divergence, the 
toroidal and poloidal fields T and S are characterized by vanishing radial component (T r = 0) and 
vanishing radial component of the curl ([V x S] r = 0), respectively 1 . 

For a self-similar problem, if there is no radial current at one r, then the same is true for all r. 
We want it to be true at large r because otherwise there would be a flow of charge to infinity, and 
the cloud would become electrically charged. Thus, in our problem, the curl of the magnetic field 
has zero radial component, and is therefore poloidal in the Stratton-Chandrasekhar terminology, 



B = S = V x 



V 



x r 



(6) 



In axial symmetry, we recover the condition B • V$ 
magnetic flux function. 



0, where <3> = — 27rsin#<9 , I'/<9# is the usual 



1 Notice that this terminology is different from the one commonly adopted in astrophysics, where is customary to 
call poloidal a field with components only along r and 9, and toroidal a field with only a component along <p. 
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3.2. Nondimensional self-similar variables 

It is easy to see from eq. (l)-(3) and eq. (6) that if p, V and ^ are chosen to have appropri- 
ate power-law dependences in r one can eliminate all r-dependences in the MS equations for an 
isothermal gas. Following Li Sz Shu (1996, hereafter LS96) we adopt the scaling 

' : R(e,<p), (7) 



2vrGr 2 



V = 2c s 2 [(l + tf )lnr + y(^)], (8) 
2r 2 r 

9 = -^F(d,<p), (9) 

where the quantity Hq in the expression of the gravitational potential eq. (8) is a dimensionless 
constant to be specified. We also define the angular parts of the V and V 2 operators, 

Vn = ^ 7 ^ + -^ 7r , = — —— S m0— + . 10 

od sin 6 dip sin 6 ov 36 sin 6 dip z 

Inserting eq. (9) into eq. (6) we obtain the expression for the magnetic field 

B = -?S-L4r + Vnn (11) 

where we have defined 



A = -V 2 n F. (12) 



Finally, the Lorentz force is given by 



i-(V x B) x B = -J-3 [(V n F ■ V n A)r - AVqA]. (13) 

Notice that a magnetic field with no radial component of its curl and proportional to 
r _1 , as assumed here, cannot be force- free. 

Inserting expressions (7)— (9) into eq. (1)— (3), we rewrite the condition of force balance in the 
radial direction as 

H R = VqF • V Q A, (14) 
and the condition of force balance in the tangential direction as 

^V n R + RVnV + AV n A = 0. (15) 
Finally, Poisson's equation becomes 

V 2 n V = R-(1 + H ). (16) 

We now show that eq. (14)-(16) generalize to the set of equations derived by LS96 for axisymmetric 
isothermal scale-free equilibria. 
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3.3. The axisymmetric case 

Assuming d/dtp = 0, and defining (f)(9) = — sin9F', we obtain 

A(9) = 



smi 



(17) 



where a prime indicate derivation with respect to 9. With this definition, eq. (14) reduces to eq. (13) 
of LS96, 



d ( 4>' 



Poisson's eq. (16) reads 



d9 V sin 9 
1 d 

sin( 



sin OdO 



dV 

_ 

d6 



zj R ■ 
-Hn— sin I 



R-l-Hn. 



(18) 



(19) 



Eliminating dV/d9 using eq. (15), and simplifying the result with the help of eq. (18), we obtain 
eq. (12) of LS96, 



1 d 



sin ( 



sine I 2H -- - 



I? 



2(R-1-H Q ), 



(20) 



that completes the set of equations governing axisymmetric, scale-free, isothermal equilibria. 
Pure hydrostatic equilibria are described by the single equation 



1 d 
ski#d# 



sin^) :2(\-R). 



(21) 



obtained by setting = and H = in eq. (20). Medvedev & Narayan (2000, hereafter MN00) 
have found an analytical solution of this equation, 



R(9) 



(liecosfl) 



2 ' 



(22) 



where iso-density contours (described by r(9) oc y 7 R{9) according to eq. [7]) are confocal ellipsoids 
of eccentricity e, with < e < 1. For e = 0, this solution reduces trivially to the singular isothermal 
sphere, with R=l. We will return to this family of equilibria in Sect. 5. 



4. Method of solution 

At first sight, the problem appears to be overconstrained, since there are four unknown func- 
tions (R, V, A and F) and four equations (eq. [12], [14], [15] and [16]), one of which (eq. [15]) has 
two components. However, it is easy to see that the vectors V^i?, and Vn^4 in eq. (15) are 

parallel, thus the condition of force balance in the tangential direction reduces to one constraint 
only. To see this, first take the curl of eq. (15), obtaining VqR x = 0, a condition implying 

that isodensity and equipotential surfaces are coincident. Then take the vector product of eq. (15) 
with V n R, obtaining VqR x VqA = 0. It follows then that V^^ x VqA = (CVD). The meaning 
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of the condition V^V x V&A = can be understood by writing the expression for the electric 
current j from Ampere's law using eq. (11), 

Eq. (23) shows that the electric current j is perpendicular to VqA, and therefore implies that the 
current flows over equipotential (or isodensity) surfaces. 

The method adopted in this work is based on the expansion of the vector variables of the 
problem in vector spherical harmonics. Any vector quantity is represented by a convergent infinite 
sequence of complex coefficients, and the problem is formulated in the Hilbert space. The solution 
of the nonlinear partial differential equations of the problem is thus reduced to the solution of 
infinite set of nonlinear algebraic equations. The coupling coefficients are expressed in terms of 
Wigner 3j symbols, and are evaluated by the Racah's formula (see Appendix A). 



4.1. Multipole expansion 

Vector spherical harmonics (see e.g. Arfken 1985) provide the natural basis for a multipole 
expansion of eq. (14)-(16). Here we follow Morse & Feshbach (1953) defining 

P/m = Y lm r, B /m = 1 VnY lm , C lm = -f x B lm . (24) 

V + 1) 

As we will see, the properties of the product of vector spherical harmonics allow non-linear terms 
to be dealt with in a systematic way by Wigner 3j symbols (or Clebsch-Gordan coefficients). 

Clearly, an infinite number of orthonormal bases can be generated from the P; m , B/ m , Q m basis 
by the application of unitary transformation. The set (24) is especially convenient for situations 
where a preferred radial direction is present, like in our case, since the P; m are radial whereas the 
B; m and C; m are tangential to the unit sphere. Thus, if we expand R, V, A and F in spherical 
harmonics, and we eliminate F and V by using eq. (12) and (16), the remaining equations of the 
problems for R and A are naturally expressed in the basis (24), with Vn^4 and VqR expanded in 
series of B; m and j in series of Q m . For the reasons explained in Sect. 4, we expect therefore the 
expansion along Q m of the equation of force balance (eq. [1]) to be trivially satisfied. 

We expand the functions R, V, A and F in spherical harmonics, 

R(e,ip) = l + H + J2RlmYim(d,v), (25) 

Im 



V{0,<p) = Y,Vi m Y lm (9,<p), (26) 

Im 

A(9, V ) = J2 A imYi m (0,<p), (27) 



Im 
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F(9, l p) = J2 F imYim(0,f), (28) 

Im 

where the sum is for / > 1 and —l<m<l. In general, Ri m , Vi m , A\ m and Fi m are complex 
coefficients. Since R, V, A and F are real functions, we have to require that 

i?,_ m = (-l) m Rt m , etc. (29) 

The constant factor l + H$ in the expansion of R(6, 92) is chosen to simplify Poisson'e equation (16). 
Notice that with this choice 

<j> Rdn = 4ir(l + H ), (30) 

a condition equivalent to the "integral constraint" of LS96 (their eq. [18]). Therefore, as in LS96, 
Ho measures the fractional increase in the mean density that arises because the magnetic field 
contributes to support the cloud against self-gravity. 

It is straightforward to compute the spherical mass-to-flux ratio, or the ratio of the mass 
contained in a sphere centered on the origin to the magnetic flux through a circle of the same 
radius in the equatorial plane. The mass enclosed in a sphere of radius r is 



o 

M - 



2vrG 



f 2c 2 r 

<b RdQ= -JL-(l + flo). (31) 



The magnetic flux through a circle of the radius r in the equatorial plane is equal by Gauss theorem 
to the magnetic flux through a semisphere of radius r, 

2c 2 r 2c 2 r r r /2 

d> = _^_ jB.rdS = -^=j)d<pJ A(e,<p)sm9 dO, (32) 
where we have used the expression of B given by eq. (11). With the expansion (27) we obtain 

47r C 2 r /-I-/ 2 yrc 2 r 

<& = -^-J>oy o ^osin0d0=^Jj>4o, (33) 

where 

* = r(i-K + iO- <34> 

Thus, in non-dimensional units, the spherical mass-to-flux ratio results 

A, = W5£=£±!» (35) 
<P ^ (piAi 

Using known properties of the vector spherical harmonics (see e.g. Varshalovich, Moskalev, & 
Khersonskii 1988), vector quantities like VqR, VqV, etc. are immediately expressed as 

V n R = Y J VWTT)Ri m Bi m , V n V = Y J VW+^)Vi m B lm , etc. (36) 

Im Im 
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In addition, the expansion in spherical harmonics presents the advantage that it makes possible 
to solve immediately Poisson's equation (eq. [16]) and the relation between A and F (eq. [12]), as 
both equations involve the angular part of the Laplacian operator. Since 

V^ m = -1(1 + l)Y lm , (37) 

eq. (12) gives the relation between the coefficients F[ m and Ai m , 

A lm = 1(1 + l)F lm , (38) 

whereas Poisson's equation (eq. 16) gives the relation between the coefficients V/ m and Ri m , 

R lm = -1(1 + l)V lm . (39) 

Inserting the expansions eq. (25)— (27) in eq. (11) and eq. (23), and eliminating the coefficients 
F[ m using eq. (38) , we obtain the expansion in vector spherical harmonics of the magnetic field and 
the electric current, 



B 



lm 



B 



lm 



and 



Using the vector relation 



i = ^^EVWTT)A lm c lm . 



(40) 
(41) 
(42) 



VxC, m = - y/l(l + l)P lm + B lm , 
r L J 

we immediately recognize in the expression (40) for B the curl of a vector potential A given by 

A = ^=y-^=C lm . (43) 



4.2. Spectral form of the equations 

We first consider the radial component of the equation of force balance, eq. (14). In terms of 
the vector spherical harmonics this equation can be written 

V^H (1 + #0)^00 + #0 Yl R Vrn'Y Vm > = Yl 

I'm' I'm' l"m" 

where we have used eq. (38) to eliminate the coefficients F\ m . Then, we multiply each term on 
both sides by Yf m and integrate over solid angle dQ, using the known orthonormality properties of 
spherical harmonics. It is convenient to introduce the coupling coefficients a^' m ", defined as 

a^ m " ee j IW • ^l"m"YC m dn, (45) 



l'(l' + 1) 1/z 
l"(l" + l) 



M'm'M"m"^l'm' 



B/// 



l"m" 



(44) 
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involving the product of three spherical harmonics (or their derivatives). The expressions for these 
coupling coefficients are given in Appendix A. We obtain from eq. (44) the condition 



4irH (l + H )5i, + H R lm = + 



l'{l' + 1) 
I" (I" + 1) 



1/2 



a^, m A Vm ,A v , m „. (46) 



I'm' l"m" 

The value of + l)a]p™' m " for I = can be easily obtained from the formulae in Appendix A. 

(-l) m ' 



/j7j i T\„ mm 



(47) 



«=0 V47T 

and using the relation (29) between coefficients of opposite m, eq. (46) for / = simplifies to 

^|^ m | 2 =47rF (l + tf ), (48) 



lm 



whereas, for I > 1, it reads 



r(r + i) 

Z"(i" + 1) 



1/2 



„mm m a a 
a U'l" A l'm' A l"m"- 



(49) 



We now expand the equation for the tangential force, eq. (15) along the two orthogonal sets 
of vectors B; m and C/ m . With the expansions given above, eq. (15) becomes 

^ - + 7= = = — Rl'm'^l'm' + E E + l)Rl'm'Vl"m"Yl'm'Rl"m" 

I'm' 2 V ' (' + 1) /' m ' F'm" 



+ E E V Z "( Z " + l)^'m'^"m''>I'm'Bi» m « = 0. (50) 
I'm' l"m" 

We take the scalar product of this equation with B? and integrate over solid angle dCl, using the 
relations eq. (39) and (38) and defining the coupling coefficients /3[Pp /m " (see Appendix A), 



omr 
Pll'l 



VW' + 1) 



Yl'm'^l"m" ■ B;* m dQ. 



(51) 



The result is 



/(/ + !) -2( 1 + ^) ^ 



2v^Tl) 



I'm' l"m" 



+ 1) 
I" (I" + 1). 



1/2 



omm m jy p 
Pll'l" til' m' til" m" 



+ E E [ /7 ( r + W'V" + l)] 1/2 ^' m "^'m^W = 0. 



(52) 



Z'm' Z"m" 



We have noticed in Sect. 4 that the tangential component of the equation of force balance has 
no component in the direction of electric-current lines. Since we see from eq. (41) that j is expressed 
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as a series containing the harmonics C/ m , we expect therefore each coefficient of the expansion of 
eq. (50) in terms of C/ m vectors to be zero. This is verified in Appendix B. 



Summarizing, the equations of the problem are: 

--Ri 



VWTT) Rlm EE 



l (l + 1) ^ mm'm" . . /_ Q s 



I" {I" + 1) 



i'm' l"m" 

from the condition of force balance in the radial direction; 

ll± ^=^ R ^ = - E E P'C + W + W* 1 "^^ 

2 VH' + 1) l'm'l"m" 

+ EE + n A™r m Rl'm'Rl"m", (54) 

I'm'V'm" L ^ 

from the condition of force balance of force in the tangential direction; and 

El^-I 2 =47r#o(l + #o), (55) 

lm 

defining the amount of support provided by the magnetic field (monopole component of the equa- 
tion of force balance in the radial direction). In this way, the representation of the physical variables 
has been transferred from function space, in terms of 6 and ip, to the infinite-dimensional Hilbert 
space, each vector component now being the coefficient of the corresponding harmonic. The pro- 
cess is analogous to transforming from the Schrodinger to the Heisenberg description in quantum 
mechanics (it is incomplete, in that the radial dependence is still represented in function space). 

The procedure for solving the equations of MS in spectral form is to select a finite set of 
coefficients by truncating the series expansion to some I = l mSLX , setting to zero all remaining 
coefficients. As a test of the method, in the next two subsections we solve the MS equations with 
m = for various values of £ max , and compare the results with the axisymmetric solutions obtained 
by MN00 and LS96 for the hydrostatic and magnetostatic case, respectively. 



5. Hydrostatic equilibria 



Setting Hq = and A\ m = for any (l,m), we obtain the set of equations governing hydrostatic 
equilibria, 



1(1 + 1) 



2 V / Z(Z + 1 



:R lr 



EE 

V I" 



I' (I' + 1) 
l"(l" + 1) 



1/2 



Pll'l" n-Vm'UV'm"- 



(56) 



Remarkably, for / = 1 both the LHS and RHS of eq. (56) are zero. To see this, observe that for 
/ = 1 the triangular relation imply that all coefficients Pifi™' m '' vanish unless /' = I", or I' = I" ±1. 
Thus, the RHS contains terms like 

amrn'm' r>2 

Pll'l' n l'm'i \P') 
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or like 



(l' + W + 2) 



1/2 



omm : m 
Pll'l'+l 



(V + + 2) 



1/2 



omm m I r> p 



(58) 



+ 1) 

Terms like (57) are zero because a and (5 coefficients are zero when l+l'+l" is odd (see Appendix A) 
terms like (58) are zero because 

V + 2 



Pu 



mm m 
"'+1 



' omm m 
'Pll'+ll' 



(59) 



as shown in Appendix A. Thus, all terms allowed by the triangular conditions in eq. (56) with 
1 = 1, are identically zero. This implies that the coefficient R\ m (or V\ m ) is undefined, an intrinsic 
degree of freedom of the problem that we refer to as the "dipole gauge". The expansion of the 
density function R contains dipole terms proportional to Yi m (9, if) with m = 0, ±1, 

R(6,(p) = l+ ^2 RimYim(0i <f) + ••• = 1 + (ci cos 9 + cos 9 cos ip + C3 cos 9 sin <p) + ... , (60) 

m=-l,0,l 

where c±, C2 and C3 are real coefficients. The last three terms represent an eccentric distortion of 
the basic spherically symmetric equilibrium (R = 1, the singular isothermal sphere) along three 
perpendicular axes, confirming the result of MN00 that the singular isothermal sphere is neutrally 
stable with respect to dipole-like density perturbations. Computing additional terms in the series 
(60), the series converges to the function R(9,ip) for each value of Ri m . Without loss of generality 
we set m = 0, equivalent to assuming that the symmetry axis of the configuration lies in the z 
direction, and for better clarity we omit the index in the expansion terms and coupling coefficients. 
For / max = 2, we have to solve for R2 as function of R\ the quadratic equation 



2R 2 = V6(3 lu Rl + (V2/%i2 + 3^221)^1^2 + V6Rl 



with solution 



R2 



14vr 



1 




20 1 



if R x < 1. 



(61) 



(62) 



The fast decrease of the coefficients with increasing Z, at least for small R±, suggests a rapid 
convergence of the series expansion. For Z max > 2, the system of truncated equations can be easily 
solved by iteration. The procedure usually converges to a solution with a relative accuracy of less 
than 10~ 5 in 4 iterations. The series solution obtained in this way corresponds to the ellipsoidal 
solutions parametrized by the ellipticity e (see Sect. 3.3), found by MN00. We show in Fig. 1 the 
solution obtained assuming Ri = 2 for Z max = 2,3 and 4. As anticipated, the terms of the series 
expansion decrease exponentially with increasing I. For comparison we also show the corresponding 
analytical solution of MN00 (eq. [22], with e = 0.466) 

In principle, the sequence of equilibria defined by the parameter Ri, with < R\ < 00, can 
bifurcate into non-axisymmetric configurations of equilibrium. A linearization of eq. (56) gives the 
condition 



1(1 + 1) 



2v^Tl) 



l{l + l)_ 



1/2 



omOm 
Pill 



+ 



1(1 + 1) 



I 1/2 



0, 



mmO 
111 



Ri = o, 



(63) 
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for the occurrence of bifurcations with angular dependence defined by harmonics with general I, m 
along the I = l,m = sequence. However, /3[^ 0m = /3[^ m0 = since 21 + 1 is odd (see Appendix 
A), and therefore no azimutally asymmetric bifurcation occurs along the sequence of hydrostatic 
equilibria generated by the R\ term. 



6. Magnetostatic equilibria with azimuthal and equatorial symmetry 

In this section we look for highly symmetrical solutions, possessing both azimuthal and equato- 
rial symmetry. We thus set m = 0, and we impose the existence of a plane of symmetry at 9 = tt/2 
(equatorial plane). The density R, the gravitational potential V, the components Bg and B v of the 
magnetic field are symmetric with respect to the equatorial plane, whereas the radial component 
of the magnetic field, B r , is anti-symmetric. This implies that the expansion of V and R contains 
multipole terms with I even, whereas the expansion of A (and F) contains multipole terms with I 
odd. Since m = 0, all expansion coefficients are real. As in the previous section, we omit the index 
m in the expansion terms and coupling coefficients. 

We truncate the infinite system of nonlinear algebraic equations to some index l max and we 
set all coefficient of the expansion equal to zero for / > ^ max . At the lowest level of approximation, 
'max = 2, we have to solve two equations for the two coefficients A\ and R2, 

^R 2 = a 2U Al (64) 

Ai = 4wH (l + H ), (65) 
where 0:211 = — 1/V1207T. The solution is 



R 2 = -2J^(1 + H ). A 1 = V47T# (l + #o), 



(66) 



The coefficients R 2 and A± determined in this way constitute the lowest-order terms of a series 
expansion for the the functions R(6) and (f)(9), 

R(9)^^(l + H )sm 2 9, (f)(9) w ^y/3H (l + H ) sin 2 9, (67) 
and the spherical mass-to-flux ratio, from eq. (35), 



This series solution converges to the solution determined by LS96 (singular isothermal toroids). As 
in the case of the hydrostatic equilibria discussed in the previous section, the convergence of the 
series solution is rapid: the approximated expression of A r , obtained with only the first term in the 
expansion of R(9) and (f)(9), is in good agreement with the numerical values calculated by LS96 as 
function of Hq, the largest discrepancy being ~ 15% for Hq » 1. 
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To obtain a more accurate series solution, for Z max > 2 we solve numerically the equations 
of MS in spectral form with Newton's method, assuming as a first guess the analytical solution 
eq. (66). The procedure converges to a solution with a relative accuracy of 10™ 8 in 4 iterations. 
In Table 1 we list the values of the coefficients A\ and Ri, and the spherical mass-to-flux ratio for 
H = 0.5 and / max = 2, 3,4, 5 and 6. The comparison with the exact numerical solution of LS96 is 
shown in Fig. 2 and 3. 

7. Three-dimensional magnet ost at ic equilibria 

Useful information about the existence of solutions of the MS equations with lower degrees 
of symmetry than those considered in Sect. 5 and 6, can be obtained by linearization of the spec- 
tral equations. The occurrence of neutrally stable configurations along the sequence of equilibria 
controlled by the parameter Hq is indicated by the vanishing of the determinant of the linearized 
system for some value of the control parameter (in our case the parameter Hq measuring the rela- 
tive amount of support provided by the magnetic field). The zeroes of the determinant signal the 
presence of neutrally stable equilibria or bifurcations points in the sense of Poincare (see e.g. Galli 
et al. 2001). 

In this way, we show in Sect. 7.1 and 7.2 that the degree of freedom represented by the 
"dipole gauge" affecting hydrostatic equilibria (Sect. 5), is also present in the case of MS equilibria 
independently on the degree of support provided by the magnetic field. The magnetic field however 
removes the degeneracy of the purely hydrostatic case, where the density distortion with I = 1 
and m = 0, ±1 gives origin to three orthogonal orientations of the same configuration. Magnetized 
equilibria instead are split by the "dipole gauge" into two families: one possessing an equatorial 
plane of symmetry but azimutally asymmetric (with m = ±1) and elongated in two orthogonal 
directions; and one without a plane of symmetry but azimuthally symmetric (with m = 0). In 
addition, the dipolar distortion of the density is coupled to a quadrupolar distortion of the vector 
potential, with the same m. 

Finally, in Sect. 7.3 we address the question of the existence of neighboring equilibria to the 
basic axisymmetric solutions determined in Sect. 6, limiting the analysis to perturbations in the 
density described by sectorial harmonics (with I = m). We find that the only allowed perturbation 
of this kind has I = m = 1, as anticipated in Sect. 7.2. 

7.1. Case I = l,m = 0: equatorially asymmetric density distortions 

First we consider axisymmetric distortions that introduce an asymmetry of the configuration 
with respect to the equatorial plane. These are represented by the small coefficients R\q for the 
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density and A 2 q for the vector potential, and are governed by the linearized equations 

V3H R 10 - V2(a%$ + 3Q?o?)i 10 ^ 2 o = 0, (69) 

and 

VSH R 10 + V2(Pffil + 3/3 1 2 1 )i? 2 oi?io - 6V2(/?? 1 ° 2 ° + f3^)A w A 2Q = 0, (70) 
where a tilde (~) indicates the coefficients of the axisymmetric solution. 

A bifurcation can occur if the determinant of the system is zero. Substituting the values of the 
coefficients, we see that the derminant vanishes for any value of Hq. For these linearized equilibria, 
the relation between i?io and A 2 q is 

H 0] J —R 10 = A 10 A 20 . (71) 

An example of these azimuthally symmetric equilibria lacking equatorial symmetry is shown in 
Fig. 4 (left panel). 



7.2. Case I = l,m = ±1: azimuthally asymmetric density distortions 

Next, we consider density perturbations azimuthally asymmetric but conserving the symmetry 
of the original equilibrium state with respect to the equatorial plane. These are characterized by 
expansion coefficients Ru and A 2 \ for the density and the vector potential, respectively, and are 
governed by the linearized equations 

VSH R U - V2(a^l + 3a 1 1 H)A 10 A 21 = 0, (72) 

and 

V3H Rn + V2(pH% + WlibRwRu - 6V2(/? 1 1 1 ° 2 1 + (3\H )Iio A 2 i = 0. (73) 
As before, the coefficients R 2 o and A\q are those of the axisymmetric solution. 

A bifurcation can occur if the determinant of the system is zero for some value of Hq. Substi- 
tuting the values of the coefficients, we see that the determinant vanishes for any value of Hq. For 
these linearized equilibria, the relation between Ru and A 2 i is 

H J^Ru = A w A 21 . (74) 

An example of these azimuthally asymmetric equilibria with equatorial symmetry is shown in Fig. 4 
(right panel). 
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7.3. Sectorial density distortions: linear analysis 

To assess the validity of Parker's theorem for axisymmetric self-gravitating equilibria, we check 
whether the class of axisymmetric solutions obtained in Sect. 6 (converging to the singular isother- 
mal toroids of LS96) allows neighbouring 3-D equilibria. To this goal, we perform a linearization 
of the set of nonlinear algebraic equations near the axisymmetric solution and we consider for sim- 
plicity distortions of the density function with arbitrary I, assuming for simplicity m = I (sectorial 
distortions). As before, the expansion coefficients corresponding to the axisymmetric solution are 
indicated by a tilde (~). 

The linearized equations for small distortions characterized by expansion coefficients Ru and 
Ai + n, with arbitrary I, are 



#0 d _ ( 2 „W /(/ + !)(/ + 2) uq 



Z(Z + l)-2(l + g ) / / 6 m 1(1 + 1) up 

2y/i{i+i) " " {yw+T) 0m+ y~e- Pm J 



_ v / 2 (/ + l)(/ + 2)(/3^ +1 + Cn)^io^+ii- (76) 
The determinant An of this systems of equations is 



2 ( l + 2 )? a l0l , am ^ , l 2 101 , /(/ + !)(/ + 2) 




A U /A 10 = H \l -^L^ +1 + /3^ n ) + I J {l + m + 2) <i + i + \l " $ ' <lu 



1(1 + 1) -2(1 + Hp) / 6 am , ll(l + l) am 



2 v / I(TTi) ly/^ + i) V 6 



Pl2l + \/ a A/2 ^20 



(77) 



and vanishes for Z = l,m = 1, as anticipated in Sect. 5 and 7.2. For I ^ 1, we have evaluated the 
determinant numerically as function of Hq for I = 2, 3, 4, 5 and 6, and < Hq < 2. The results 
are shown in Fig. 5. At least in this part of parameter space, the determinant is always positive, 
its value monotonically increasing with Hq and I, and shows no sign of having a zero for particular 
values of Hq. We can then safely conclude that the axisymmetric MS equilibria determined by 
LS96 and discussed in Sect. 6, have no other neighboring equilibria than those allowed for all values 
of Hq by the "dipole gauge" discussed in Sect. 6, 7.1 and 7.2. 

In this sense, these non-symmetric equilibria represent the only exceptions (within the assump- 
tion of the present study) to a generalized version of Parker's theorem, originally formulated for 
systems with translational symmetry, extended to self-gravitating equilibria with axial symmetry. 
In the next section we consider the equilibria originated by the "dipole gauge" from a physical 
point of view, and we conclude that all these solutions are probably not force-free at the origin, a 
singular point for a scale-free configuration, and therefore they cannot represent realistic models of 
equilibrium of isolated cosmic bodies. 
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8. Summary 

The results obtained in this paper show that previously known axisymmetric solutions of the 
MS equations for an isothermal self-gravitating gas, under the hypothesis of scale invariance and 
global neutrality, allow only neighbouring equilibria characterized by a / = l,m = 0, ±1 angular 
dependence of the density distortion for any value of the degree of magnetic support (including 
zero and infinite). For m = 0, the original axisymmetric equilibrium is distorted by a "bending" 
of the isodensity contours with respect to the equatorial plane, preserving azimuthal symmetry; 
for m = ±1, the equilibrium is distorted by a "stretching" of the isodensity contours along one 
of two orthogonal directions in the equatorial plane, preserving up-down reflection symmetry. In 
the absence of magnetic fields (for Hq = 0) these two classes of distorted equilibria reduce to the 
ellipsoidal equilibria found by MN00. In the limit of vanishing thermal support (Hq — ► oo) the 
configuration of equilibrium reduces to a thin disk supported only by magnetic tension against its 
self-gravity, and the m = 1 density distortion corresponds to the elliptical disklike equilibria found 
by Galli et al. (2001). What is the significance of the neutral stability of these configurations to 
density distortions characterized by a dipolar angular dependence? 

According to MN00, the ellipsoidal hydrostatic equilibria are not force-free at the origin, where 
the gravity of all the matter of the configuration produces a non- vanishing force trying to restore 
symmetry with respect to the equatorial plane. According to Cai & Shu (2004), the same conclusion 
holds for the azimutally asymmetric solutions found by Galli et al. (2001) for MS equilibria in the 
thin-disk limit. These two classes of equilibria stem from the singular isothermal sphere and the 
singular isothermal disk, respectively, that in turn are the limiting cases (for Hq = and Hq = oo) 
of the family of singular isothermal toroids of LS96. It is therefore tempting to conclude, that 
all singular isothermal equilibria, characterized by a distortion of the isodensity (or equipotential) 
surfaces proportional to the I = 1 harmonics, are not force-free at the origin for any value of Hq. 
If this is the case, their relevance to represent realistic equilibria is doubtful. 

The same conclusion probably hold for configurations rotating with spatially uniform veloc- 
ity Uu, (the only rotation law compatible with isothermality and spatial self-similarity). Galli et 
al. (2001) found that rotating singular isothermal disks are neutrally stable to m = 1 perturbations 
for any value of the rotation rate. In a similar vein, MN00 found that the the azimuthally and equa- 
torially symmetric rotating models of Toomre (1982) and Hayashi et al. (1982) can be "continued", 
for any value of the rotation velocity, into a sequence of axisymmetric equilibria lacking equatorial 
symmetry (in our language, originated by a I = l,m = density distortion). At the light of the 
results described in Sect 7.1 and 7.2, it is natural to expect that the sequence of Toomre-Hayashi 
models also possess non-axisymmetric counterparts with a dominant I = l,m = 1 asymmetry, 
although solutions of this kind are not known. Likely, also these hypothetical rotating asymmetric 
equilibria are not force- free at the origin. The I = 1, m = 0, ±1 distortion, the allowed distortion of 
axisymmetric equilibria found in the present study, may then represent a gauge freedom that creeps 
somehow into general self-similar isothermal equilibria, irrespectively of the presence of magnetic 
fields, or rotation, or else. The appearance of this gauge freedom in self-similar isothermal systems, 
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and its physical significance, deserves further scrutiny. As a counterexample, it should be easy 
to show that, assuming a non-isothermal (e.g. polytropic) equation of state, the analogues of the 
Li-Shu magnetized equilibria studied by Galli et al. (1999), (or the polytropic analogues of the 
Toomre-Hayashi models) are not affected by this gauge freedom. 

For arbitrary I = m 7^ 1, a perturbation analysis shows that the symmetric magnetized solu- 
tions found by LS96 do not have neighbouring non-axisymmetric equilibria. This results is analo- 
gous to Parker's theorem for systems with translational invariance, and suggests that the validity 
of Parker's theorem can be extended to self-gravitating axisymmetric equilibria. Combining this 
result with the findings of Galli et al. (2001), one is led to the conclusion that the presence of 
fluid motions (specifically, super-magnetosonic rotation) is a crucial ingredient for the occurrence 
of symmetry break-ups in MHD equilibria. 

The results of this paper imply that, under the assumed conditions, very few (and possibly 
not physically meaningful) non-axisymmetric solutions of the steady MHD equations do exist. 
Of the assumed conditions, probably the most severe are the assumptions that (i) the magnetic 
field is analytic everywhere, and that (ii) the new, non-symmetric solutions are accessible from 
the basic states by regular perturbations (i.e. small-parameter expansions). As for the former 
assumption, Moffatt (1985) has shown that magnetostatic equilibria of non-trivial topology in a 
perfectly conducting fluid may generally contain tangential discontinuities (current sheets), i.e. 
they cannot be described in terms of analytic functions; in contrast, equilibrium fields that are 
analytic functions of space are subject to severe structural constraints, as shown by Arnol'd (1965, 
1966). As for the latter assumption, our results show that equilibria of non-trivial topology, if they 
exist, cannot in general be reached starting from axisymmetric states by regular perturbation. This 
is consistent with the conclusions of Rosner & Knobloch (1982), that the response of stationary 
solutions of nonlinear equations to finite-amplitude, symmetry-breaking perturbations may not in 
general be obtained in terms of small-parameter expansions of the variables. These limitations 
should be kept in mind in interpreting the results of this paper. 



8.1. Implications for molecular cloud cores 

Recent statistical studies (Jones, Basu, & Dubinski 2001; Jones & Basu 2002; Goodwin, Ward- 
Thompson, & Whitworth 2003) based on available catalogues of molecular cloud cores and Bok 
globules (typical size L = 0.1 pc, sound speed c s = 0.2 km s _1 , average density n(H2) = 10 4 cm~ 3 , 
and typical magnetic field strength B = 10 /^G), do not support the possibility that cores are 
axisymmetric configurations. A good fit to the observed axial distribution is generally found as- 
suming instead that cores are triaxial ellipsoids. The best-fit axial ratios a : b : c determined 
statistically for molecular cloud cores (1 : 0.9 ± 0.1 : 0.5 ± 0.1 according to Jones & Basu 2002; 
1 : 0.8 ± 0.2 : 0.4 ± 0.2 according to Goodwin, Ward-Thompson, & Whitworth 2003), suggest that 
cores are preferentially flattened in one direction and nearly oblate (a « b > c), and imply that 
they may not be particularly far from conditions of equilibrium. 
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Taking these observational results at face value, one is led to consider the fate of a cosmic 
cloud with its frozen-in magnetic field formed by whatever process in a configuration lacking a high 
degree of symmetry and presumably not in an exact equilibrium state. As discussed in Sect. 2.2, 
Parker (1979) argued that realistic magnetic fields with no well defined symmetries must evolve in 
a genuinely time-dependent way, until all non-symmetric components of the field are destroyed by 
dissipation and reconnection and the topology becomes symmetric. This process can hardly be of 
any relevance for the interstellar gas, where ohmic dissipation times are larger than the age of the 
Universe. Ambipolar diffusion on the other hand can only redistribute the mass inside flux tubes 
and drive a stable equilibrium to the threshold of dynamical instability, but cannot dissipate the 
magnetic energy stored in the field. 

To the extent that interstellar clouds can be represented as isolated MS equilibria, these sys- 
tems must instead undergo Alfven oscillations (weakly damped by ambipolar diffusion) around the 
equilibrium state with period (Woltjer 1962) 

(78) 



(c 2 + 4)V 2 ' 

where L is the size of the system, c s is the sound speed and v\ is the Alfven speed. This kind 
of behaviour is evident in the numerical and analytical calculations of Hennebelle (2003) rela- 
tive to the homologous evolution of prolate and oblate magnetized isothermal spheroids. In the 
non-selfgravitating case, the stability of magnetostatic equilibria of arbitrary complex topology 
was studied by Moffatt (1986) through construction of the second variations of the magnetic and 
kinetic energies with respect to a virtual displacement about the equilibrium configuration. Mof- 
fatt's (1986) results show that a general class of space-periodic magnetostatic equilibria is stable to 
disturbances of arbitrary lengthscale. If perturbed in some way, the fluid responds executing oscil- 
lations about this equilibrium, that are eventually damped if due account is taken of viscosity. If 
we adopt the deviation from axial symmetry of the shapes of cloud cores as rough measure of their 
nonequilibrium, a simple harmonic oscillator analogy provides an estimate of the average velocity 
(u 2 ) 1 / 2 of pulsation, 

(c 2 +v 2 A ) a 2 

that, for a : b = 1 : 0.8 gives (v 2 ) 1 / 2 m 0.2 {c 2 + v 2 A ) 1 / 2 . For typical conditions of molecular 
cloud cores, this implies coherent pulsations with period t«4x 10 5 yr, average velocity (v 2 ) 1 ^ 2 « 
0.05 km s _1 , and maximum velocity ~ 0.1 km s _1 . These pulsational motions may have already 
been detected in the isolated globule B68 (Lada et al. 2003). 



It is a pleasure and a privilege to thank Frank Shu, who proposed the problem, suggested the 
solution technique, and anticipated the results, long before the calculations presented in this paper 
were completed. An anonymous referee is also aknowledged for a very insightful report that helped 
to improve the original manuscript. 
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A. The coupling coefficients 

The coupling integrals of eq. (45) and (51) can be expressed in terms of integrals of three 
spherical harmonics using standard recurrence formulae to eliminate the derivatives. The resulting 
integrals can be evaluated with the Gaunt formula, 

jYC m Y Vm ,Y vlm „ dn = [1(1 + 1)1' (I' + 1)1" (I" + lp 2 N wl „Cfflp; m ", (Al) 



where 



and 



Nun" = - 



(2l + l)(2l' + l)(2l" + 1) 
4vrZ(Z + l)l'(V + 1)1" (/" + !) 



1/2 



(A2) 



^"-(-ir(o o o)(- m I !!»)■ (A3) 

The 3j symbols are algebraically defined by e.g. the Racah formula (see e.g. Landau & Lifshitz 1997, 
Varshalovich, Moskalev & Khersonskii 1998). The evaluation of the coupling coefficients is best 
achieved by group-theoretical methods and the use of Wigner 6j and 3j symbols. Here we quote 
only the results, details on the procedure can be found in Jones (1985). In terms of the quantities 
Nun" and G^' m ", we have 

„mm'm" a at /~rmm'm" nmm'm" a at /~imm! m" I K A\ 

a ll'l" = A l^ll'l"^U'l" , Pll'l" = A l>Nll>l"<^U>l>> , ( A4 ) 

where 

A, = I" (I" + 1) + I' (I' + 1) - 1(1 + 1), k v = I" (I" + 1) + 1(1 + 1) - I' (I' + 1). (A5) 

The coupling coefficients a^5p' m " and (3lfi™' m " are real, and are equal to zero when the following 
(triangular) conditions on the Wigner 3j symbols are not all satisfied: 

l<l' + 1", I' <l + 1", I" <l + 1', l + + 1" even, (A6) 

and 

m = m! + m". (A7) 

We recall that a 3j symbol is invariant under even permutation of columns and is multiplied by 
the phase factor (— l)( l+l ' +l ") for odd permutations. Thus G^,™,'™ 1 " being the product of two 3j 
symbols, is invariant for any permutation. Using this property, one obtains, for example, 



7 _j_ o 

nmm'm" omm" m' i a q\ 

Plll+l - — —Pll+ll ■ l A «V 
The coefficients o^™' m " and Pul™' m " are real, whereas 7[p™' m " is imaginary: 

„mm'm" --p at Tjmm'm" I 

lll'l" ~ ~ lV l^WV'Hivi" ) ( Ay ) 
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where 
and 

= [(/' + 1" - + 1 - /")(-/' + 1 + 1" + + i' + i + 1)] 1/2 . (An) 

Notice that all the dependence of the coupling coefficients on the azimuthal number m is contained 
in the Wigner 3j symbols, whereas the remaining factors depend on I only (the Wigner-Eckart 
theorem). 



B. Proof that one component of the equation of force balance is trivially satisfied 



To show this, we take then the scalar product of eq. (50) with C* lm and integrate over solid 
angle, obtaining the condition 



EE 



+ 1) i 1/2 



^ ^ i" a" + 1) 

I'm' l"m" 1 v ' A 

where we have defined 



TIP/" m Rl'm'Rl"m" ~ E E + + l)] 1 ^ 2 7lr™ m M'm'M"m" = 0, 

(Bl) 



I'm' l"m" 



_ ,111-111/ lie 

lll'l" = 



(B2) 



The expressions of the coupling coefficients 7[p™ m in terms of the Wigner 3j symbols is given in 
Appendix A. It is not difficult to see that, as expected, this equation is trivially satisfied, since the 
two terms in eq. (Bl) are equal to the coefficients of the expansions of VqR x VqV and VqA x VqA 
in terms of P/ m , both of which are zero. 

We begin by expanding the vector product VqR x VqV, 

V n R x V n V = Y,Y1 Ri'm'Vi» m "[l\l' + 1)1" \l" + l)] 1/2 (B rm ' x B,» m »). (B3) 

I'm' l"m" 

Eliminating the coefficients Vi» m " with eq. (39), and using the relation (see Jones 1985), 

B// m / x Bj» m « = -(B// m / • Q» m »)r, (B4) 

we obtain 



V n R x V n V = J2Y1 



'm' l"m" 



/'(*' + 1) 
I" (I" + 1) 



1/2 



Rl'm'Rl"m"(Bl 



(B5) 



We now compute the expansion of this equation in series of Pi m . For this, we need to evaluate the 
coefficients 

^m B 2'm' • Cl"m" dVL = <h [Yi m B*, m , ■ C*„ m „r dtt 
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= (-1)™' j[Y lm B v _ ml ■ c?„ m „r dn = (-ir' v ^7Ti)[7Sr i '- m > 

= VW+^)bwn" ml T = -VW+^)bwu' m "T = VW+^hu^'" 1 ", 



(B6) 



where we have used the definition of 7^™ m and the symmetry properties of Wigner 3j symbols. 
Thus, we finally obtain 



VnR x V n V = £ + l)P, m i £ £ 



+ 1) 



1/2 



(I'm' l"m" 1 v ; J 

The same procedure, applied to the cross product VqA x VqA, shows that 



JW™ m Rl'm'Rl"m" ( ■ (B7) 



V^A x V n A = - £ 7^Tl)P, m I £ E [*'(*' + W'V" + l)] 1/2 7j?r m "^'m'^i"m" [ • (B8) 



lm 



.I'm' l"m" 



Since VqR x V^V = 0, as can be shown by taking the curl of eq. (15), and VqA x VqA = 0, 
we conclude that for each (I, m) the quantities inside curly brakets in eq. (B7) and (B8) must be 
zero. Thus, each term of the expansion of the equation of force balance along Q m , eq. (Bl), is 
zero. 



Table 1. Multipole Coefficients for H = and m = 



'mm -Rl R-2 R3 Ra 



2 2 0.7847 

3 2 0.7631 0.2604 

4 2 0.7623 0.2546 0.0792 
exact 2 0.7623 0.2542 0.0791 



Table 2. Multipole Coefficients for H q = 1/2 and m = 



^max 




R2 


A 3 


i? 4 


A 5 


Rq 




2 


3.0670 


-2.3780 










2.002 


3 


3.0583 


-3.2861 


-0.2670 








1.943 


4 


3.0591 


-3.2573 


-0.2581 


0.6839 






1.944 


5 


3.0593 


-3.2496 


-0.2532 


0.8388 


0.0356 




1.940 


6 


3.0593 


-3.2496 


-0.2534 


0.8249 


0.0324 


-0.1523 


1.941 



Fig. 1. — Iso-density contours for the axisymmetric hydrostatic (Hq = 0) equilibrium with R\ = 2 
obtained with the method described in this paper {long- dashed curves, Z max = 2; short-dashed 
curves, / max = 3; dotted curves Z max = 4). The analytical solution of MN00 with the same value of 
R\ is shown by the solid curves. 
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Fig. 2. — Density function R{9) and flux function (f)(6) for the axisymmetric case with Hq = 0.5. 
The solid curves are the numerical solutions of LS96. The solutions obtained with the method 
described in this paper for m = are shown by long-dashed curves (Z max = 2), short-dashed curves 
(^max = 4), and dotted curves (/ max = 6). The lower panels show the differences between the exact 
solution of LS96 and the series solution obtained in this work. 



Fig. 3. — Axisymmetric magnetostatic equilibrium for H$ = 0.5 (isodensity contours and magnetic 
field lines). The solution of LS96 is shown by solid lines, the solution obtained in this paper with 
^max = 6 is shown by dotted lines. The agreement of the two solutions is very good. 
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1=1 m=0 1=1 m = ±l 




Fig. 4. — The I = 1 gauge. For any value of Hq, magnetostatic equilibria with azimuthal and 
equatorial symmetry (dashed lines) possess neighboring equilibria with a non-zero I = 1 density 
component (solid lines). The two panels show examples of these equilibria with small density 
perturbations containing I = l,m = and I = l,m = ±1 harmonics. 
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Fig. 5. — The determinant Au of the system of linearized equations describing sectorial (I = m) 
distortions of axisymmetric magnetostatic equilibria. The determinant is normalized to A\q, the 
value of the dipole term in the expansion of the vector potential for the magnetic field. Notice that 
the determinant is zero for any value of Hq for I = l,m = ±1 density distortions, as discussed in 
Sect. 7.2, and non-zero in all other cases. 
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